Rapid inverse radiative transfer solver for multiparameter spectrophotometry without integrating sphere

Abstract. Significance Multiparameter spectrophotometry (MPS) provides a powerful tool for accurate characterization of turbid materials in applications such as analysis of material compositions, assay of biological tissues for clinical diagnosis and food safety monitoring. Aim This work is aimed at development and validation of a rapid inverse solver based on a particle swarm optimization (PSO) algorithm to retrieve the radiative transfer (RT) parameters of absorption coefficient, scattering coefficient and anisotropy factor of a turbid sample. Approach Monte Carlo (MC) simulations were performed to obtain calculated signals for comparison to the measured ones of diffuse reflectance, diffuse transmittance and forward transmittance. An objective function has been derived and combined with the PSO algorithm to iterate MC simulations for MPS. Results We have shown that the objective function can significantly reduce the variance in calculated signals by local averaging of an inverse squared error sum function between measured and calculated signals in RT parameter space. For validation of the new objective function for PSO based inverse solver, the RT parameters of 20% Intralipid solutions have been determined from 520 to 1000 nm which took about 2.7 minutes on average to complete signal measurement and inverse calculation per wavelength. Conclusion The rapid solver enables MPS to be translated into easy-to-use and cost-effective instruments without integrating sphere for material characterization by separating and revealing compositional profiles at the molecular and particulate scales.


Introduction
Conventional spectrophotometers are probably the most commonly used type of instruments for material analysis by determination of absorbance A or attenuation coefficient μ t as a function of wavelength λ for a given sample.For turbid materials including biological tissues, additional optical parameters are needed to characterize molecular composition and particle sizes on scales close to λ by their ability to absorb and scatter light.One paradigm is to retrieve these parameters by the radiative transfer (RT) theory, 1 and its spectroscopic implementation is termed multiparameter spectrophotometry (MPS).The RT theory quantifies light-matter interaction by the parameters of absorption coefficient μ a and scattering coefficient μ s in addition to a single-scattering phase function pðθ; ϕÞ with θ and ϕ as the polar and azimuthal scattering angles.MPS thus requires measurement of multiple light scattering signals to solve multiple inverse scattering problems (ISPs) at selected values of λ that remains challenging for highly turbid samples.2][3] Alternative model functions have been studied for pðθ; ϕÞ but p HG ðcos θÞ is the preferred one with a form fully specified by an anisotropy factor g as the mean value of cos θ.By choosing p HG as the phase function, MPS is to measure multiple light scattering signals and determine inversely the RT parameters as function of at λ.These parameters can be expressed as a vector of PðλÞ ¼ ðμ a ; μ s ; gÞ.
][5][6][7][8][9][10][11] If only one integrating sphere is used, signals need to be measured in two steps with one for R dh and another one for T dh with T c acquired in either step.To solve the ISPs, one first derives μ t (¼ μ a þ μ s ) from T c by the Beer-Lambert law followed with retrieval of μ s and g from R dh and T dh using different methods of forward modeling guided by gradient descent based inverse solver.The modeling methods include numerically solving the RT boundary-value problem, adding-doubling algorithm or Monte Carlo (MC) simulations.Despite its popularity as a research tool, the need for two integrating spheres or two steps of signal measurement makes it difficult to translate this approach into an easy-to-use instrument like a conventional spectrophotometer.In addition, the use of integrating sphere limits severely the accessibility of the approach to non-specialists due to time-consuming sample assembly and system maintenance.Other approaches without integrating sphere have been investigated such as goniometric measurement and detection of T c and non-hemispherical R d and T d at fixed angles with MC based inverse algorithms. 12,13e have previously shown that the RT parameter vector P can be uniquely determined from three simultaneously measured signals of non-hemispherical diffuse reflectance R d , diffuse transmittance T d and forward transmittance T f without integrating sphere. 14,15MC simulations are performed to accurately calculate signals as R dc ; T dc , and T fc and repeated to match the measured ones.A gradient descent algorithm has been developed to guide iteration in the RT parameter space to minimize an objective function δðPÞ defined as the sum of squared percentage differences between the measured and calculated signals.An ISP at λ is deemed as solved with P s when δðP s Þ ≤ δ th in which δ th represents a threshold based on the experimental errors in signal measurement.The inherent statistical variance in the calculated signals, however, makes it difficult to accurately solve ISPs for highly turbid and optically thick samples by gradient decent despite the existence of a unique solution. 16Specifically, sizable regions of small δ values exist in the RT parameter space of P for these ISPs of large scattering albedo að¼ μ s ∕μ t Þ and optical thickness τð¼ μ t DÞ with D as sample thickness.In such regions, values of δðPÞ fluctuate considerably due to the variance of MC simulations that often lead to errors in the solution given by P s .To solve these challenging ISPs for highly turbid samples, one needs to either significantly increase the number of photons in MC simulations for variance reduction or search manually by contour analysis of δðPÞ distributions in the RT parameter space.Either way gives rise to high computational cost and prevents rapidly solving ISPs for MPS. 16In this report, we present a rapid inverse solver for MPS based on a particle swarm optimization (PSO) algorithm with a novel objective function ρ 1∕δ ðPÞ through local averaging in the space of P to determine P s ðλÞ. 17The function ρ 1∕δ ðPÞ significantly reduces the effect of MC simulation variance in calculated signals on inverse calculation by PSO and computation time to solve ISPs by reducing the number of photons in MC simulations.Our validation results with 20% intralipid samples demonstrate that P s ðλÞ can be retrieved rapidly between 520 and 1000 nm by executing MC simulations on one GPU board.

Signal Measurement by MPS
An experimental system has been constructed to measure R d ðλÞ, T d ðλÞ, and T f ðλÞ for validation of the new inverse solver with the signal detection configuration shown in Fig. 1(a).The details of the experimental system were previously reported. 18Briefly, a xenon light source (XL1-175-A, WavMed Technologies Corp.) and a monochromator (CM110, CVI Corp.) are employed to produce a monochromatic beam with λ adjustable between 520 and 1000 nm in steps of 20 nm and bandwidths around 5 nm.The beam is modulated at a frequency of f 0 ¼ 370 Hz by a mechanical chopper (SR540, Stanford Research Systems) and incident on an assembly consisting of a turbid sample confined in a spacer ring between two glass slides.The intensity of the incident beam I 0 is monitored by a photodiode of D 1 (FDS1010, Thorlabs, Inc.).Three photodiodes (FDS100, Thorlabs, Inc.) of D 2 to D 4 are used to measure respectively I Rd for diffusely reflected, I Td for diffusely transmitted and I Tf for forwardly transmitted light intensity.The current signals of photodiodes were amplified by an in-house built four-channel lock-in amplifier to obtain the measured signals of R d ¼ I Rd ∕I 0 ; T d ¼ I Td ∕I 0 and T f ¼ I Tf ∕I 0 .Figure 1(a) shows the detection configuration for acquisition of measured signals from the sample assembly.

Signal Calculation by iMC
An in-house developed individual photon tracking MC (iMC) code was employed to calculate signals as R dc ; T dc and T fc from given P for a phantom of the same shape as the sample inside a spacer between two glass slides by tracking N 0 photons, which imports the parameters of sample size and detection configuration as input data. 16,19,20The code injects each of the N 0 photons incident on the sample assembly and then tracks the photon once it transports inside a glass slide or sample until it is either absorbed inside the sample or escapes into air.The exit location and propagation direction of an escaping photon on a glass slide surface are used to determine if it hits a detector for detection.A counter associated with each detector records the number of detected photons as N Rd ; N Td , and N Tf by the detector D 2 ; D 3 , and D 4 , respectively.The above process repeats until the total number of injected photons reaches N 0 and the calculated signals are given by R dc ¼ N Rd ∕N 0 , T dc ¼ N Rd ∕N 0 and T fc ¼ N Tf ∕N 0 .Because of independence in trajectory among the N 0 photons, the numbers of detected photons follow Poisson distributions. 21To estimate the variance in these photon numbers, one can draw a random number q of Poisson distribution with fðq; q m Þ as probability mass function and q m as the mean.In the case of N Rd calculated by an iMC simulation, q m equals to N 0 R dc∞ if distribution of N Rd follows fðN Rd ; q m Þ and R dc∞ yields the variance-free value of R dc by tracking "infinite" number of photons.One thus can use fðq; q m Þ to quantify the effect of variance on calculated signals and optimize the objective function for variance reduction.

PSO Algorithm
A stochastic algorithm based on PSO has been developed to guide iMC simulations and solve for P s ðλÞ in the RT parameter space from measured signals by optimizing an objective function as shown in Fig. 1(b).An objective function quantifies the difference between measured signals and calculated ones obtained by given PðλÞ, and optimization of this function reduces the difference to solve for P s ðλÞ from an initial choice of PðλÞ.The PSO algorithm was chosen for its high efficiency and ability to perform global search or avoid local traps in the RT parameter space.A search proceeds through multiple threads in PSO, and the threads are represented by a swarm of "particles" with index i ∈ ½1; I and I as the number of threads or particles.Here, a particle i symbolizes a thread of positions in the RT parameter space along which it moves from Pði; jÞ to Pði; j þ 1Þ as 22 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 4 ; 6 0 4 where j ∈ ½1; J; J is the number of iterations; v is the particle's velocity; q and q' are random numbers of uniform distribution from 0 to 1; P p ðiÞ is the particle-best position; P g is the swarmbest position; w 0 ; w p , and w g are weights of terms contributing to vði; j þ 1Þ; and χ is a damping constant to increase stability.A particle transits toward P p ðiÞ and P g through completed iterations by Eq. ( 1) that respectively yields the best position of that particle and all particles for optimized objective function to reduce difference between calculated and measured signals.By setting values of w p and w g equal and large in comparison to w 0 in Eq. ( 1), for example, enables fast convergence of all particles from their best positions P p ðiÞ towards P g on average as j increases.A region Γ s in the RT parameter space is designated as the search space and Pði; j þ 1Þ in Eq. ( 1) is set to Pði; jÞ if the former moves out of Γ s .The search stops once the objective function at P g reaches a preset threshold or j exceeds J and current P g is saved as P s ðλÞ for output.
3 Results and Discussion where δ f ðPÞ ¼ jS c∞ ðPÞ − 1j 2 is the variance-free value of δðPÞ.It is clear that δ f ðP s Þ ¼ 0 if P s represents a "true" solution of ISP.

Comparison of Two Objective Functions for Variance Suppression
Based on above analysis, local averaging on δðPÞ in the RT parameter space was first employed to suppress variance and hence increase the stability of inverse solutions.We then compared two objective functions defined below on their abilities to further reduce the effect of variance E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 7 ; 6 5 1 where A P is a set of m P vectors consisting of P and its neighbors on a 3D grid of P. The components of f c as defined in Eq. ( 2) can be obtained as ðq∕q m Þ − 1 by drawing a random number q independently from fðq; q m Þ and setting the q m values, related to R dc∞ ; T dc∞ , and T fc∞ , equal to focus on the effect of variance for simplicity.We further set m P ¼ 7 with P and its nearest neighbors as the members of A P and assumed that all of them have the same value of δ f ðPÞ.The values of ρ δ ðPÞ and ρ 1∕δ ðPÞ are plotted in Fig. 2(a) against δ f ðPÞ ranging from 0.1% to 1% with a solid line representing δ f ðPÞ or δ f ðPÞ −1 .These results show clearly that ρ 1∕δ ðPÞ fluctuates significantly less than ρ δ ðPÞ for δ f ðPÞ between 0.6% and 1%, which is the typical range of δ th set by the error level of measured signals. 16o gain further insight on fluctuation in calculated signals, the relative difference Δ between ρ δ ðPÞ and δ f ðPÞ or ρ 1∕δ ðPÞ and δ f ðPÞ −1 are calculated and averaged as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 4 7 4 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 7 ; 4 1 9 where δ f is uniformly sampled between 0.1% and 1.1% in each sum and m δ is the total number of δ f samples in the sums.The sampled values of δ f ðPÞ correspond to the mean values of each signal term in δ f ðPÞ ranging from 1.8% to 6.1% since δ f ðPÞ is defined as the sum of squared relative differences of calculated and measured values on three signals.Figure 2(b) illustrates the dependence of the two Δ functions on q m values with m δ set to 1000.As discussed before, q m represents the mean value of fðq; q m Þ which is the probability mass function of Poisson distribution for modeling detected photon numbers, such as N Rd and related calculated signal R dc : Therefore, the ratio of q m ∕N 0 equals to R dc∞ as the variance-free value of R dc in the above case and larger q m corresponds to larger N 0 .The results in Fig. 2(b) demonstrate that use of ρ 1∕δ ðPÞ instead of ρ δ ðPÞ or δðPÞ as the objective function can significantly reduce the variance of calculated signals for the same value of N 0 .Alternatively, N 0 can be considerably reduced for the same variance to speed up forward calculations if ρ 1∕δ ðPÞ is adopted.For example, computational time of iMC simulations can be cut in half with ρ 1∕δ ðP), instead of ρ δ ðPÞ, as the objective function by tracking half of photons with similar variance.

Measurement of 20% Intralipid Sample
To validate the new approach, RT parameters of 20% intralipid (I141-100ML, Sigma-Aldrich) have been determined from three measured signals of R d ; T d , and T f for λ from 520 to 1000 nm in steps of 20 nm.S2 in the Supplementary Material. 23The same measurements were repeated on another sample of D ¼ 81 μm and the inversely determined RT parameters agree well with those in Fig. 3 within the experimental errors.
The search region Γ s for PSO based inverse algorithm at each wavelength was set between 1.00 × 10 −4 and 2.00 mm −1 for μ a , 10.0 and 200 mm −1 for μ s , 0.10 and 1.00 for g.In additional to choosing I ¼ 27 for the particle number and J ¼ 100, we set w 0 ¼ 0.900, w p ¼ w g ¼ 2.05, and χ ¼ 0.730 for the PSO parameters defined in Eq. ( 1) for optimized performance of inverse calculations.The initial positions of Pði; 1Þ for all particles were randomly distributed in each of nine equal partitions of Γ s to ensure globally optimized solutions of ISPs.We found that the convergence toward the final solution was not affected by the choice of Pði; 1Þ under the above condition.After calculation of ρ 1∕δ with Pði; 1Þ, iterations as defined in Eq. ( 1) were performed in Γ s .Specifically, P p ðiÞ was set as Pði; 1Þ for particle i and P g as P p ði 0 Þ that has the largest ρ 1∕δ value among all particles for j ¼ 1 to obtain vði; 2Þ and Pði; 2Þ by Eq. ( 1).The search then continued as j increases to 2, 3,. . .and stopped when either ρ 1∕δ ðP g Þ exceeds ρ th for j ≤ J or j > J with P g saved as P s ðλÞ and output together with ρ 1∕δ ðP g Þ and δðP g Þ.The threshold ρ th corresponds to the value of ρ 1∕δ in Eq. ( 5) with δ th set to 0.8% and m P to 7.
It took about eight iterations on average to solve an ISP of P s ðλÞ from the measured signals of 20% intralipid sample.The total numbers of iMC simulations were found to be about 800 for solving each ISP if 27 particles were employed in PSO based search.Signal were calculated by iMC simulations in which N 0 was gradually increased from 5 × 10 5 to 5 × 10 6 as δ became 1% or less.On average, over different values of RT parameters and N 0 , each iMC simulation took about 0.2 s and solving P s ðλÞ per wavelength took about 160 s by executing the code on one GPU board (Nvidia, GeForce RTX 2080 Ti).The results of P s ðλÞ in Fig. 3 compare similarly to those obtained by methods using integrating sphere. 8,24For example, both of the range and wavelength dependence of μ s and g agree well with those in Ref. 8 for the 20% intralipid sample while the range of μ a is between those of Ref. 24 for 10% intralipid and Ref. 8. Since the 20% intralipid is high turbid with scattering albedo a ¼ μ s ∕μ t > 98% over the measured wavelengths, the large differences in the values of μ a can be attributed to the very large errors in determining μ a with values smaller than μ s by three orders of magnitude as previously discussed. 8We note further that the values of optical thickness τ and albedo a of the intralipid sample reach the maximum values of 7.1 and 99.98%, respectively for λ ¼ 520 nm.Results of our study on samples of 20% intralipid with sufficiently large D values showed that current MPS method with the inverse solver reported here can yield unstable solutions of P s when the value of τ becomes larger than 7.1 for very large values of a above 99.9990%.A detailed comparison of the measured and calculated signals and analysis of background noise in signal measurement suggests that the instability of inverse solution is mainly caused by the variance in calculated signals, which may be mitigated without using N 0 much larger than those used for ISPs of τ ≤ 7.1 for iMC simulations.A study is underway to further improve the PSO based inverse solver by taking into account the δðPÞ distribution obtained from completed iterations, which is expected to enable accurate measurement of RT parameters for highly turbid samples with τ up to 15 without significant increase of computational cost.Such an improvement can make the approach of MPS concerned here capable of characterizing optically thick samples by their RT parameters for which the conventional approach with integrating sphere fails for inability to accurately measure collimated transmittance T c .

Conclusions
We have analyzed the effect of variance by MC simulations on solving ISPs and proposed a novel objective function to reduce variance in calculated signals and speed up simulations.Combination of the objective function with the PSO algorithm leads to a rapid inverse solver for determination of RT parameters of turbid samples from three measured signals without integrating sphere.The inverse solver has been validated by the results of RT parameter retrieval on samples of 20% intralipid in a wavelength range of 520 to 1000 nm with ISPs solved rapidly using one GPU board.Taken together, we have demonstrated a new approach of MPS, which allows its translation into a powerful and easy-to-use instrument for clear separation and characterization of molecular composition and turbidity.

Disclosures
Authors declare no conflicts of interest.

Fig. 1
Fig. 1 (a) Configuration of signal detection with incident beam indicated by the red line: D: sample thickness; D 1 : for monitoring I 0 ; M: mirror; D 2 to D 4 : for measurement of I Rd ; I Td , and I Tf ; θ 0 : angle of incident beam from z-axis; d R ; d T , and d f : distance between the origin and front center of D 2 , D 3 , and D 4 ; and θ R ; θ T , and θ f : angle of sensor surface normal (blue dash lines) of D 2 , D 3 , and D 4 from the z-axis.(b) Work-flow chart of signal measurement in brown boxes and inverse calculation in green boxes.

Figure 3
presents the RT parameters of a sample with thickness D ¼ 102 μm obtained by the PSO based algorithm with ρ 1∕δ ðPÞ selected as the objective function.Signal measurement was repeated three times to obtain the mean values and standard deviations which are plotted in Fig. S1 in the Supplementary Material.Previously determined values of real refractive index n r ðλÞ of the 20% intralipid were used to obtain interpolated values in iMC simulations at the wavelengths of measurement for the calculated signals and n r ðλÞ is presented in Fig.

Fig. 3
Fig. 3 The wavelength dependence of RT parameters of one 20% intralipid sample of D ¼ 102 μm in thickness.The error bars were inversely determined from different combinations of the mean and standard deviation values of the measured signals shown in Fig. S1 in the Supplementary Material on five selected wavelengths.The red line is a power law fitting of μ s ¼ Cλ −2.340 with C ¼ 1.626 × 10 8 .
where R dc∞ , T dc∞ , and T fc∞ are the variance-free signals calculated by tracking "infinite" numbers of photons.It should be noted that f c is independent of P and one can express S c ðPÞ by f c 3.1 Effect of Variance in Calculated Signals on Solving ISPsMPS requires solving one ISP for each value of λ from the measured signals.The first step is to perform forward calculations of signals from given P ¼ ðμ a ; μ s ; gÞ by iMC simulations of lightmatter interaction in the sample.The calculated signals can be normalized by respective measured signals and expressed as a vector of S c ðPÞ ¼ ðR dc ðPÞ∕R d ; T dc ðPÞ∕T d ; T fc ðPÞ∕T f Þ.An inverse algorithm is to iterate iMC simulations toward the objective of making S c ðPÞ as close to 1 ¼ ð1;1; 1Þ as possible.A function of squared error sum δðPÞ ¼ jS c ðPÞ − 1j 2 is often selected as the objective function to guide iteration and decide when stops.Ideally, δðP s Þ vanishes for a perfect ISP solution given by P s .In practice, one deems an ISP as solved if δðPÞ ≤ δ th with δ th representing a threshold determined by the error level in signal measurement.To quantify the effect of variance in calculated signals by a stochastic iMC simulation, we define a fluctuation vector for calculated signals as with S c∞ ðPÞ ¼ ðR dc∞ ðPÞ∕R d ; T dc∞ ðPÞ∕T d ; T fc∞ ðPÞ∕T f Þ and "o" denoting the vector form of component-wise product.Combining these results, δðPÞ can now be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 4 ; 1 2 4